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Abstract. Combining density-matrix and Lanczos algorithms we propose a new 
optimized phonon approach for finite-cluster diagonalizations of interacting electron- 
phonon systems. To illustrate the efficiency and reliability of our method, we inves- 
tigate the problem of bipolaron band formation in the extended Holstein Hubbard 
model. 



1 Introduction 

Considerable work is currently focused on the experimental and theoreti- 
cal study of strongly coupled electron-phonon (EP) systems, triggered by 
the recognition that the EP interaction plays an important role in under- 
standing the physics of novel materials such as colossal magneto-resistance 
manganites |Q or the very recently discovered superconducting magnesium 
diboride @. From a theoretical point of view the challenge is to describe the 
partly exotic properties of these materials in terms of simplified microscopic 
models which take into account the complex interplay of charge, spin and 
lattice degrees of freedom. 

As a generic model for systems with competing electron-electron and 
electron-phonon interactions the extended Holstein Hubbard model (EHHM), 

(i.j):<7 i i,l;cr i 

(1) 

is usually considered ||,fi|5|], where c^} and b^ annihilates [creates] a spin-u 
electron and a phonon at Wannier site i, respectively, and n ir7 — c irr c ilT . The 
Hamiltonian ([[]) consists of a kinetic term describing the electronic motion 
on a discrete lattice (transfer amplitudes t), an extremely screened (on-site) 
Coulomb repulsion (Hubbard parameter U), and a "density-displacement" 
type non-screened EP coupling (oc kxq) 



= m _ j|2 + 1)3/2 ' X ° = VWMUO , KXQ = y/SpUJQ . (2) 
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Here ujq is the bare phonon frequency of dispersionsless optical phonons, 
being polarized in the direction perpendicular to the chain (ID case). Defining 
the polaron binding energy as e p = (xq/ojq)J2i fi : (0) — l-27e p , the famous 
Holstcin Hubbard model (HHM) results by setting 

fi(i) = K,6 ti i , e p —> e p , (3) 

i.e., with respect to the EP coupling term the EHHM represents an extension 
of the Frohlich model [|| to a discrete ionic lattice or of the Holstein model 
including longer ranged EP interactions. 

Adapting the EHHM to real physical situations one is frequently faced 
with the difficulty that the energy scales of electrons (t, U), phonons (u>o) 
and their interaction (e p ) are of the same order of magnitude, causing an- 
alytic methods, and especially adiabatic techniques, to fail in most of these 
cases. Thus, at present, the most reliable results came from powerful numeri- 
cal calculations, such as finite-cluster exact diagonalizations (ED) ||,|9|,[l0| or 
(Quantum) Monte Carlo simulations which are usually performed on 

supercomputers. But even for these numerical approaches strong EP inter- 
actions are a demanding task, since they require some cut-off in the phonon 
Hilbert space. Starting with the work of White Q in 1993, during the last 
years a class of algorithms became very popular, which based on the use of 
a so-called density matrix for the reduction of large Hilbert spaces to man- 
ageable dimensions. 

In the present paper, we will demonstrate that finite-cluster diagonaliza- 
tion methods also benefit substantially from these ideas. Along this line, in 
the next section we introduce an optimized phonon approach for the ED of 
electron-phonon problems. To exemplify this technique, we analyze the for- 
mation of bipolarons in Sec. III. Our main results are summarized in Sec. IV. 

2 Optimized phonon approach 

Let us first resume the connection between density matrices and optimized 
basis states. Starting with an arbitrary normalized quantum state 

D..-1 D v -1 

M= E E ^\v)\r) (4) 

r=0 !/=0 

expressed in terms of the basis {|^)|r)} of the direct product space H = 
H u <g) H ri we wish to reduce the dimension D v of the space H v by introducing 
a new basis, 

Dv-l 

!/=0 



Density-Matrix Algorithm for Phonon Hilbert Space Reduction 



3 



with v — . . . (Dy — 1) and D„ < D v . The projection of \ip) onto the corre- 
sponding subspace H = Hp ® H r C H is given by 



D r -1 D D -\ D„-l 



1$ = Y Y Y °^u'1 v -t\^)Y) 

r=0 0=0 v'=0 
D r -1 Dp-1 D v -1 

= Y Y Y ^vv^lv'lv'MV) ■ (6) 
r— />— v,v'= 

We call {p)} an optimized basis, if \ip) is as close as possible to the original 
state Therefore we minimize |||V ; )~lV')l| 2 with respect to the elements aa v 
of the transformation matrix a under the orthogonality condition (v'\v) = 
X^=o 1 °*v'v a vv = &v'v- Applying the latter condition, we find 

\\Wi - i^)ii 2 = w#> - <v#> - (m + «#> 

D r -1 D„-l D„-l 

Y Y Y ^r^«L'7„'r+ Ec - 
r— v— v.v'— 

£> 5 -l D u -1 

7'— O^'—O v.v' ,v" — 
_D r -l D^-l D^-l 

r— z>— i/.z^'— 

= l-T^apa*), (7) 

where p — 1 itrlfr is called the density matrix of the state |^) with 

respect to {|^)}- We observe immediately that the states {\v)} are optimal if 
the rows of a are eigenvectors of p corresponding to its Dc, largest eigenvalues 

Following Zhang et al. |Q, we now apply these features to construct an 
optimized phonon basis for the eigenstates of an interacting electron/spin- 
phonon system. Consider a system composed of N sites, each contributing 
a phonon degree of freedom \vi), Vi = O...00, and some other (spin or 
electronic) states \ri). Hence, the Hilbert space of the model under consid- 
eration is spanned by the basis {^^Sq 1 1^)1^)}- Of course, to numerically 
diagonalize an Hamiltonian operating on this space, we need to restrict our- 
selves to a finite-dimensional subspace. To calculate, for instance, the lowest 
eigenstates of the HHM (|l])-(||) we could limit the phonon space spanned by 
= (^iO _1/ ' 2 (^l) I/ *|0) by allowing only the states Vi < Di. Most simply 
we can choose Di = M V i yielding D p h = M for the dimension of the 
total phonon space. However, if we think of the states {(Q^Jq 1 |z^)} as eigen- 
states of the Hamiltonian Tt p h — ^Eilo 1 ^!! ^ ^ s m ore suitable for most 
problems to choose an energy cut-off instead. Thus we used the condition 
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J^^q 1 Vi < M, leading to D p h = ( N+ ^^ 1 ), for most of our previous numer- 
ical work (see e.g. Ref. JTj|). For weakly interacting systems already a small 
number M of phonon states is sufficient to reach very good convergence for 
ground states and low-lying excitations. However, with increasing coupling 
strength most systems require a large number of the above 'bare' phonons, 
thus exceeding capacities of even large supercomputers. In some cases one 
can avoid these problems by choosing an appropriate unitary transformation 
of the Hamiltonian (e.g. by using center-of-mass coordinates), but in general 
it is desirable to find an optimized basis automatically. 



large site 




small sites □ 

□ □ □ □ □ 



1 ~ i ■■■■ ■■■■ ■■■■ ■■■■ ■■■■ 

Fig. 1. Structure of the phonon basis in terms of the highest accessible /ij. 
Left: as proposed by Zhang et al. [|l4j; Right: used within this work. 

Within the current density-matrix algorithm fll4| for the construction of 
an optimal phonon basis the phonon subsystem is considered as a product 
of one 'large' (i = 0) and a number of 'small' sites (i > 0). Each site except 
the large one uses the same optimized basis {|Mi>o)} = {|^)} with v = 
. . . (m — 1), while the basis of the large site consists of the states {\v)} plus 
some bare states {|f)}, {l^o)} = ON({|i>)} U {|^)}), where ON(. . . ) denotes 
orthonormalization (see Figure [l]). After a first initialization the optimized 
states are improved iteratively through the following steps 

(1) calculating the requested eigenstate \ip) of the Hamiltonian H in terms 
of the actual basis, 

(2) replacing with the most important (i.e., largest eigenvalues wo) 
eigenstates of the density matrix p, calculated with respect to and 

{|M0», 

(3) changing the additional states {\v)} in the set {|/xo)}> 

(4) orthonormalizing the set {|/^o)}; arid returning to step (1). 

A simple way to proceed in step (3) is to sweep the bare states {\v)} through 
a sufficiently large part of the infinite dimensional phonon Hilbcrt space. 
One can think of the algorithm as 'feeding' the optimized states with bare 
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phonons, thus allowing the optimized states to become increasingly perfect 
linear combinations of bare phonon states. Of course the whole procedure con- 
verges only for eigenstates of TL at the lower edge of the spectrum, since usu- 
ally the spectrum of a Hamiltonian involving phonons has no upper bound. 
The applicability of the algorithm was demonstrated in Ref. [Q with the 
Holstein model (i.e., U — in Eq. (jl])) as an example. 

When we implemented the above algorithm together with a Lanczos ED 
method for our systems of interest, we found two objections against the above 
choice of an optimized basis: (i) the basis is not symmetric under the sym- 
metry operations of the Hamiltonian (e.g. translations), and (ii) the phonon 
Hilbert space is still large (-D p h = M m , where M is the dimension at the 
large site), since we usually need more than one optimized state per site. 

The first problem is solved by including all those states into the phonon 
basis that can be created by symmetry operations (see Figure [lj right panel) , 
and by calculating the density matrix in a symmetric way, i.e., by adding 
the density matrices generated with respect to every site, not just site i = 0. 
Concerning the second problem we note that the eigenvalues wa of the density 
matrix p decrease approximately exponentially, see Figure 0. If we interpret 
Wi, ~ exp(— av) as the probability of the system to occupy the corresponding 
optimized state \v), we immediately find that the probability for the complete 
phonon basis state (^^g 1 is proportional to exp(— a^^Q 1 This is 
reminiscent of the energy cut-off discussed above, and we therefore propose 
the following choice of phonon basis states at each site, 

Vi: {\IH)} = ON({| M )}) (8) 
I ^ | opt. state \v), < fj, < m 

1 bare state \v), m < /i < M ' 

and for the complete phonon basis |®x , i /i i <M I J"*) } > yielding D p h = ( N+ x~ X ) ■ 

To discuss the nature of the obtained optimized states, the convergence 
of the algorithm, and some variants in more detail, let us consider a special 
case of the HHM, Eq. (0), namely the Holstein model of spinless fermions in 
one dimension, 

n = -* E [ c l c i+i + H - c ] + 9^0 5>! + - 1) + E b t b i ■ ( 10 ) 

i i i 

The optimized phonon approach comes into play in the case of strong electron- 
phonon coupling g. Then the systems develops lattice distortions which ac- 
company the itinerant fermions. These finite elongations need to be expressed 
in terms of Harmonic oscillator states \v t ) = {vJ)~ 1,2 {b\) Vi |0), that are cen- 
tered around the equilibrium position. Hence, for the numerical diagonal- 
ization either a large number of these 'bare' states or some other states em- 
bodying a finite distortion are required. By 'sweeping' through the large space 
of 'bare' phonons, the above optimization procedure automatically creates a 
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Fig. 2. Convergence of the ground-state energy for the Holstein model, 
Eq. ([n]), with electron-phonon coupling g — 5 and different phonon fre- 
quencies: (a) luo — 0.1 1, (b) ujq = t, and (c) lu = 10 1. Solid lines: one 
local set, dotted lines: two local sets for each fermion state, dot-dashed lines: 
momentum dependent sets. 



small number of basis states \ which are sufficient for a good approximation 
of eigenstates of H(t, g, ujq)- 

In Figure ^| the convergence of the ground state energy of a two-site system 
at half-filling is shown for an increasing number of iterations (solid lines). 
We compare the results of an ordinary diagonalization using up to M = 80 
bare phonons per site (i.e., Z? p h = ( N+ ^~ 1 ) — 3240) and of the optimized 
approach. In the latter case the phonon basis consists of 6 optimized and 4 
bare states, i.e., M = 10 and _D p h = 55. Each optimized state is chosen to be 
a linear combination of the first 120 bare states. Initially we set \v) — \ v) and 
then sweep the 4 bare states through the states \v) with v = . . . 119. Vertical 
dashed lines denote the end of each sweep. The plateau structure is due to 
the fact that states of high v are less important for the optimized basis \v). 
Note that every iteration involves the calculation and diagonalization of the 
density matrix, an update of the operators op , which need to be transformed 
to the current basis, and, most expensive, a Lanczos iteration to obtain a new 
approximation for the requested eigenstate of Tt. It is therefore recommended 
to use the optimized states obtained for a small cluster as the initial basis 
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Fig. 3. Eigenvalues wp of p calculated with the ground state of the Holstein 
model for g = 5 and frequencies: (a) ujo — 0.1 1, (b) w = t, and (c) uj = Wt. 
Filled circles: one local set; open squares: two local sets for each fermion state; 
filled diamonds and triangles: momentum dependent sets. 



for a larger cluster, and to restart the Lanczos procedure with the previous 
eigenvector (although it is expressed in a slightly different basis). 

The figure also includes data for two variants of the algorithm, namely the 
construction of two sets of optimized states, one for each local fermion state, 
or a transfer of the calculation into momentum space and the use of different 
optimized states for each momentum q. The advantages of these ideas become 
clear looking at Figure ^|. Here the eigenvalues wp of the density matrix p 
are given for different phonon frequencies loo and coupling g = 5. In the anti- 
adiabatic regime of high frequencies we observe pairs of equally important 
eigenstates of p [filled circles in panels (b) and (c)]. This indicates that the 
lattice follows the fermion immediately (small polaron) . It is locally distorted 
to one side or the other, if the site is occupied by an electron or not. Hence, 
only half of the optimized states couple to one of the local fermion states, 
and it appears reasonable to use different basis sets for the two situations. 
This results in step free exponential decrease of the eigenvalues of the density 
matrices for both basis sets (open squares in Figure which allows for a 
smaller cut-off M in the total phonon space of the cluster. 

However, at small frequencies the lattice is slow compared to the fermions, 
and distortions are long ranged. Obviously, there is no gain using different 
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local basis sets, the convergence as well as the decay of the eigenvalues Wo of 
the density matrix is slow (see panel (a) in Figures || and |J) . The performance 
of each local optimization seems to be poor and switching to momentum 
space is recommended. After a shift of operators, b+ — ► &| + § , and Fourier 
transform the Hamiltonian 7i, Eq. (fio|), reads 

n = - i E [ c k+i + H - c -] + E( & i + + w ° E 6 i & i + ^ ( n ) 

= -2* E cos [f fc ] + 5= E( fet -* + + w ° E fot A + ^ > 

k k,q k 

with £J S = g 2 uJo (J2k nk ~ t) ■ Using the same cut-off in the total phonon ba- 
sis as described before, we now optimize each g-mode separately. The results 
turn out to be excellent for all phonon frequencies. The ground-state energy 
converges quickly (dot-dashed lines in Figure |J) and only a few optimized 
states are required, as can be seen from the rapidly decreasing eigenvalues of 
the density matrix (filled diamonds and triangles in Figure [5j) . 

For illustration, in Figure ^ we give optimal wave functions obtained for 
the frequencies loq = 0.1 1 and 10 1 with v increasing from top to bottom. 
In the left panels bold lines denote the wave functions of a single local set, 
whereas solid and dashed lines mark the two sets of functions that depend on 
the fermion occupation number. The right panels show optimal wave func- 
tions in momentum space. In all cases x denotes the normalized elongation 
-^{b^ + b), i.e., the expansion of the optimized states \u) — ~Y^ U Oi~ v \v) is 
plotted using Harmonic oscillator eigenfunctions in elongation space, 

-x 2 /2 

H„(x), (12) 
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with Hermite polynomials H v (x). In momentum space only one or two states 
have a non- negligible eigenvalue wo, therefore higher optimized states, which 
do not contribute to the ground state of H, are not expected to be con- 
verged. In all cases the most important states (y = 0) resemble the eigen- 
states of shifted Harmonic oscillators. If we use a single local set, the states 
correspond to symmetric and anti-symmetric combinations of left- and right- 
shifted oscillator functions (remember the step structure in Figure |J). For 
ojq = 0.1 1 the shift of the real space functions decreases with higher v. This 
reflects the fact, that the lattice is slow and fermions move within a distor- 
tion. In comparison, the case of high phonon frequency luq = 10 1 is much 
simpler, as we deal only with states of almost fixed shift. 

In summary, the simple example of the two-site Holstein model provides 
good insight into the properties of the optimized phonon approach. It is 
very efficient for determining an optimal basis within a given phonon Hilbert 
space. Nevertheless, we are not relieved from choosing the most appropri- 
ate decomposition of our model under consideration. Here physical intuition 
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and some knowledge of the model is necessary to decide between real and 
momentum space, or other special choices for the phonon Hilbert space. As 
demonstrated above, the structure of the optimized wave functions and the 
behaviour of the eigenvalues of the density matrix may give some hints. 



3 Bipolaron band formation in the EHHM 

Besides bi- /polaron formation itself, the question of whether polarons or bipo- 
larons can move itinerantly has been the subject of much controversy over 
the last decades (see, e.g., the debate in Ref. |1(|). The existence of polaronic 
bands has been verified by ED techniques in ID and 2D |L7||l8|Jl0| , however, 
since both the width and the (electronic) spectral weight of the polaronic 



bands are exponentially reduced in the strong-coupling case 19 , the coher- 
ent band motion of small bi-/polarons becomes rapidly destroyed, e.g. by 
impurities or thermal fluctuations. 

Recently it has been discovered that a longer-range EP interaction leads 
to a decrease in the effective mass of polarons ||,|| and bipolarons || in 
the strong-coupling regime, which can have significant consequences because 
the quasi-particles are more likely to remain mobile. Indeed, for the single- 
electron extended Holstein model (Q) the polaron band dispersion was shown 
to be less renormalized as compared to the Holstein model |I| . Here we would 
like to present some results for the two-electron case, where the competition 
between attractive EP interaction and Coulomb repulsion becomes impor- 
tant. 

By applying the optimized phonon approach outlined in Sec. 2 to EHHM 
we obtain the lowest eigenvalues of the two-particle system in each K sec- 
tor (E2(K)), presented in Fig. The gain in performance compared to or- 
dinary ED is illustrated in Table |l|. The EP parameters e p /t = 3.0 and 
LOo/t = 0.5 are chosen in order to address the intermediate coupling and 
frequency regime, which is almost impossible to investigate analytically. At 





phonon 
cut-off 


matrix 
dimension 


Lanczos 
diaganilizations 


memory 
requirement 


CPU time 
per run 


ED 


M = 81 


~5x 10 13 


1 


> 1000 TBytes 


??? 


optimized 
phonons 


M = 13 
m = 7 


1.8 x 10 6 


50-300 


~ 500 MBytes 


~ 20 - 200 
CRAY T3E hrs 



Table 1. Problem sizes and computer requirements for the exact diago- 
nalization method (ED) and the optimized phonon approach to solve one 
parameter set of the bipolaron problem on a ten-site lattice within the same 
accurancy. For both strategies the calculation of the ground-state in a fixed 
K-sector is assumed. 
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first we note that the two electrons always form a bipolaronic bound state in 
the EHHM: The bipolaron binding energy A — E 2 {0) — 2£i(0) is negative, 
irrespective of the Hubbard interaction strength U (see inset of Fig. ^). This 
is an important difference between the EHHM and the HHM. In the HHM, 
a critical interaction U c exists for any EP coupling where the bipolaron un- 
binds Then, calculating the if-dependent binding energy, defined as 
A(K) = E 2 (K)-mm K ,, K ,,[Ei(K')+Ei(K")] with K' + K" = K{mod 2tt), 
we find a bound bipolaron for all if-values. This means that the dispersion 
curves depicted in Fig. |^ can be deemed to be well-defined bipolaron quasipar- 
ticle bands. Of course, due to the "dressing" with phonons, the effective mass 
of the bipolaron is substantially enhanced. Accordingly the coherent band- 
width of the EHHM bipolaron, AE — £2(71") — £2(0), is reduced as compared 
to the free electron case but, on the other hand, it is notably larger than those 
of the HHM bipolaron, indicating that the EHHM bipolaron is a rather mo- 
bile quasiparticle. As U increases AE monotonously decreases, which clearly 
is a correlation effect (put in mind that U hinders double occupancy). At 
last we notice that the band structure of the bipolaron significantly deviates 
from a rescaled cosine band only for U much smaller than 2e p . In this case 
the on-site and nearest-neighbour density-density correlation functions are 
about the same size. For U > 2e p the band becomes almost cosine shaped. 
This striking cosine band dispersion, indicating free-particle like behaviour, 
was previously observed for the HHM model near U — 2e p (Refs. [ p^plj ], cf. 
also open symbols in Fig. ^) , and has been attributed to the formation of an 
inter-site bipolaron. If the inter-site bipolaron moves as a bound pair through 
the lattice, owing to the retardation effect (ujq < t), the second electron can 
take the advantage of the lattice distortion left by the first one, still avoid- 
ing the on-site Coulomb repulsion. As a result the residual bipolaron-phonon 
interaction is small. 

To gain more insight into the differences between EHHM and HHM bipo- 
laronic states, we have calculated the bipolaron kinetic energy, E\a n , given by 
the ground-state average of the first term of (|j) . Figure || presents £kin as a 
function of the EP interaction parameter e v at fixed U/t = 6 and wo/t = 0.5. 
Comparing the results for the EHHM and HHM, we can distinguish three 
regimes. For small EP couplings, the two polarons are bound in the EHHM, 
but not in the HHM. The higher kinetic energy of the EHHM bipolaron is 
mainly due to the larger incoherent contribution to the f-sum rule, which orig- 
inates from incoherent hopping processes of the two charge carriers within 
the joint lattice distortion spread over the whole lattice. In the intermedi- 
ate coupling regime predominantly inter-site bipolarons are formed in both 
models. Now the coherent part to the f-sum rule being proportional to the 
(inverse) effective mass) plays a decisive role. The EHHM bipolaron has to 
drag a larger phonon cloud than the HHM bipolaron coherently through the 
lattice and therefore acquires a larger effective mass. At strong EP couplings, 
in the HHM, a second transition to an on-site bipolaronic state takes place 
at about e p /t = 3(~ U/2), whereas the EHHM bipolaron stays in a spatially 



Density-Matrix Algorithm for Phonon Hilbert Space Reduction 11 

much more extended state. Consequently we observe a very gradual decrease 
of the bipolaron kinetic energy for the EHHM. Finally we comment on the 
renormalization of the coherent bandwidth of the EHHM bipolaron. As in 
the single polaron case, at small EP couplings the band structure is flattened 
near the zone boundary by the intersection with the dispersionsless phonon 
branch. Hence we find AE ~ 0.5 for e p /t <C 1. The strong reduction of the 
bandwidth starting to come in at about e p /t > 1 can be traced back to the 
formation of a bipolaron with pronounced nearest-neighbour correlations. 

4 Summary 

The objective of this work was the presentation of an advanced phonon opti- 
mization algorithm for application in Lanczos diagonalizations. At the heart 
of our procedure an 'optimized' phonon basis is created automatically im- 
proving the most important eigenstates of the density matrix by means of 
gradual 'sweeps' through a sufficiently large part of the infinite dimensional 
phonon Hilbert space. Within this scheme the density matrix is calculated in 
compliance with the symmetries of the underlying model. Depending on the 
physical problem under consideration, the efficiency of the proposed method 
might be improved considerably using more than one set of optimized phonon 
states or by performing the optimization procedure in momentum space. 

The reliability of our approach was demonstrated by calculating the band 
dispersion of bipolarons in the framework of the extended Holstein Hubbard 
model. In comparison with the Holstein Hubbard model where with increas- 
ing strength of the EP interaction a sequence of transitions from two unbound 
large polarons to an inter-site bipolaron and finally to a self-trapped on-site 
bipolaron takes place, in the EHHM the two electrons form a bipolaronic 
bound state for all EP couplings, irrespective of the magnitude of the Hub- 
bard interaction. As an effect of the longer ranged non-screened EP coupling 
included in the EHHM, the EHHM bipolaron is a rather mobile spatially 
extended quasiparticle even in the extreme strong-coupling regime. 
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Fig. 4. Optimized phonon states in elongation space for (a) loq = 0.1 1, and 
(b) ujq = lOt, and different choices of the optimized basis (left panels: real 
space; right panels: momentum space). 
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Fig. 5. Bipolaron band structure in the ID extended Holstcin Hubbard 
model. Data points (filled symbols) are obtained from exact diagonalizations 
of the EHHM on eight- and ten-site chains (employing periodic boundary 
conditions) with different Hubbard interaction strengths; dotted lines give 
the corresponding rescaled cosine bands having the same bandwidths. At 
U/t = 6, the bipolaron band dispersion of the Holstcin Hubbard model is 
included for comparison (open symbols). The inset displays the bipolaron 
binding energy A as a function of the inverse Hubbard interaction strength. 
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Fig. 6. Rcnormalizcd bipolaron kinetic energy in the eight-site EHHM (filled 
symbols) and HHM (open symbols). The inset gives the coherent bandwidth 
of the bipolaron in dependence on the EP coupling strength. 



